{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Surface restructuring analysis of LAMMPS trajectory - NEB preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Header"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lmpdump as lmpdump    # Use lmpdump.py parser\n",
    "import os\n",
    "import math\n",
    "import numpy as np\n",
    "import time\n",
    "from tqdm import tqdm\n",
    "\n",
    "print('Please cite DOI: 10.1021/acs.jpcc.9b04863.')\n",
    "print('This script detects interlayer surface restructuring events from LAMMPS MD trajectory and generates images for subsequent LAMMPS NEB calculation of each event.')\n",
    "print('This script is interactive and requires user input. Please read carefully before proceeding:\\n')\n",
    "\n",
    "print('Note 1: Intended only for interlayer surface restructuring events in FCC metals.')\n",
    "print('Note 2: Intended only for periodic slab models with vacuum along the z-direction.')\n",
    "print('Note 3: Intended only for LAMMPS NVT simulations.')\n",
    "print('Note 4: Must have used DFT unit cell if planning on subsequent DFT optimization (via NEB2DFT.py).\\n')\n",
    "\n",
    "print('Note 5: Requires a LAMMPS DATA file containing equilibrated box dimensions.')\n",
    "print('Note 6: Requires two trajectories, ordered by atom ID & unscaled: (1) Dumped every ps, xy-unwrapped; (2) Dumped every 0.1 ps, wrapped.')\n",
    "print('Note 7: You may want to clamp the first trajectory to ideal FCC lattice sites (via lmpclamp.py) before event detection.\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load LAMMPS trajectories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwd = os.getcwd()\n",
    "\n",
    "# LAMMPS custom-style dump file\n",
    "name = input('Enter the name of the LAMMPS trajectory file (ID, type, xu, yu, z) dumped every ps : \\nWarning: Must be ordered by atom ID, unscaled & xy-unwrapped! Please quit now if not so.\\n')\n",
    "print('Loading trajectory-1...')\n",
    "file = pwd + '/' + name\n",
    "xyz1 = lmpdump.lmpdump(file, loadmode='all')\n",
    "print('\\n')\n",
    "\n",
    "print('Loading time steps...')\n",
    "step = []\n",
    "for s in xyz1.finaldict.keys():    # keys = time steps\n",
    "    step.append(s)\n",
    "print('Done.\\n')\n",
    "\n",
    "size = np.array(step).size    # Number of time steps\n",
    "N_dump1 = xyz1.finaldict[0][1].shape[0]    # Number of atoms dumped\n",
    "xlo_ext = xyz1.finaldict[0][0][0]    # x_ext = x + xy\n",
    "xhi_ext = xyz1.finaldict[0][0][1]\n",
    "\n",
    "# LAMMPS atom-style dump file\n",
    "name = input('\\n Enter the name of the LAMMPS trajectory file (ID, type, x, y, z) dumped every 0.1 ps : \\nWarning: Must be ordered by atom ID, unscaled & wrapped! Please quit now if not so.\\n')\n",
    "print('Loading trajectory-2...')\n",
    "file = pwd + '/' + name\n",
    "xyz2 = lmpdump.lmpdump(file, loadmode='all')\n",
    "N_dump2 = xyz1.finaldict[0][1].shape[0]    # Number of atoms dumped\n",
    "print('\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract unit cell information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load LAMMPS data file after NPT equilibration\n",
    "name = input('Enter the name of the LAMMPS DATA file containing the equilibrated box dimensions : ')\n",
    "file_eq = pwd + '/' + name\n",
    "log = open(file_eq,'r')\n",
    "log_lines = log.readlines()\n",
    "log_lines = [line.split() for line in log_lines]\n",
    "\n",
    "for line_index, line in enumerate(log_lines):\n",
    "    if line:\n",
    "        for l in line:\n",
    "            \n",
    "            if l == 'atoms':\n",
    "                N_atom = int(line[0])\n",
    "                \n",
    "            elif l == 'types':\n",
    "                N_type = int(line[0])\n",
    "\n",
    "            elif l == 'xlo':\n",
    "                xlo = float(line[0])\n",
    "                xhi = float(line[1])\n",
    "            \n",
    "            elif l == 'ylo':\n",
    "                ylo = float(line[0])\n",
    "                yhi = float(line[1])\n",
    "            \n",
    "            elif l == 'zlo':\n",
    "                zlo = float(line[0])\n",
    "                zhi = float(line[1])\n",
    "            \n",
    "            elif l == 'xy':\n",
    "                xy = float(line[0])\n",
    "                xz = float(line[1])\n",
    "                yz = float(line[2])\n",
    "            \n",
    "            elif l == 'Atoms':\n",
    "                head = line_index + 1\n",
    "                start = line_index + 2\n",
    "            \n",
    "            elif l == 'Velocities':\n",
    "                end = line_index - 1\n",
    "\n",
    "Lx = xhi - xlo    # Box dimensions\n",
    "Ly = yhi - ylo\n",
    "Lz = zhi - zlo\n",
    "tan = Ly/xy\n",
    "\n",
    "xyz = []    # N*3 matrix\n",
    "for line in log_lines[start:end]:\n",
    "\n",
    "    line[:2] = np.array(line[:2]).astype(int)    # Atom ID, atom type\n",
    "    line[2:5] = np.array(line[2:5]).astype(float)    # x, y, z\n",
    "    line[5:8] = np.array(line[5:8]).astype(int)    # ix, iy, iz (inverse sign)\n",
    "\n",
    "    x = line[2] + line[5]*Lx + line[6]*xy    # Wrapped coordinates\n",
    "    y = line[3] + line[6]*Ly\n",
    "    z = line[4] + line[7]*Lz\n",
    "    \n",
    "    xyz.append([x,y,z])\n",
    "xyz = np.array(xyz)\n",
    "\n",
    "ID2 = int(input('Enter the ID of a nearest-neighbor atom that lies in the same layer as atom #1 : '))\n",
    "ID3 = int(input('Enter the ID of another nearest-neighbor atom that lies in the same layer as atom #1 : '))\n",
    "ID4 = int(input('Enter the ID of an atom that lies one layer above atom #1 : '))\n",
    "\n",
    "for atom, coord in enumerate(xyz):\n",
    "    \n",
    "    if atom+1 == 1:\n",
    "        atom1 = coord    # atom1 = Reference atom\n",
    "    \n",
    "    elif atom+1 == ID2:\n",
    "        atom2 = coord    # atom2 = 1st coplanar atom\n",
    "\n",
    "    elif atom+1 == ID3:\n",
    "        atom3 = coord    # atom3 = 2nd coplanar atom\n",
    "        \n",
    "    elif atom+1 == ID4:\n",
    "        atom4 = coord    # atom4 = Atom one layer above\n",
    "\n",
    "r12 = np.linalg.norm(np.subtract(atom1,atom2))\n",
    "r13 = np.linalg.norm(np.subtract(atom1,atom3))\n",
    "dr_avg = np.mean([r12, r13])    # In-plane NN distance\n",
    "\n",
    "a1 = np.subtract(atom2,atom1)    # 1st intralayer vector\n",
    "a1 = a1 / np.linalg.norm(a1)\n",
    "\n",
    "a2 = np.subtract(atom3,atom1)    # 2nd intralayer vector\n",
    "a2 = a2 / np.linalg.norm(a2)\n",
    "\n",
    "a3 = np.cross(a1,a2)    # Normal vector\n",
    "if a3[2] < 0:\n",
    "    a3 = -a3    # Ensure normal vector along +z\n",
    "a3 = a3 / np.linalg.norm(a3)\n",
    "\n",
    "dn_avg = np.dot(atom4,a3) - np.dot(atom1,a3)    # Interlayer distance; n = normal component\n",
    "\n",
    "log.close()\n",
    "\n",
    "print('Unit cell information & normal direction extracted.\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Event detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Parsing trajectory for restructuring events...')\n",
    "t_begin = time.perf_counter()\n",
    "\n",
    "file_ev = pwd + '/event_flag.txt'    # List of time steps flagged as event points\n",
    "out_ev = open(file_ev,'w')\n",
    "\n",
    "header1 = ' tf = final time (flagged step) \\n dt = event duration \\n ni = initial normal component \\n nf = final normal \\n CNi = initial coordination number \\n CNf = final coordination number \\n \\n'\n",
    "header2 = 'tf (ps) \\t dt (ps) \\t ID \\t type \\t ni \\t nf \\t CNi \\t CNf \\n'\n",
    "out_ev.write(header1)\n",
    "out_ev.write(header2)\n",
    "\n",
    "skip_index = [0, 1, 2, size-2, size-1]    # Skip time steps at the edges of the trajectory\n",
    "\n",
    "for s_index, s in enumerate(step):\n",
    "    \n",
    "    if s_index in skip_index:\n",
    "        continue\n",
    "\n",
    "    ev_atom = []    # List of active atoms detected\n",
    "    for atom in range(0,N_dump1):\n",
    "            \n",
    "        # f0 = current time step (final 0)\n",
    "        # f1 = 1 step after (final 1)\n",
    "        # f2 = 2 steps after (final 2)\n",
    "        # i1 = 1 step before (initial 1)\n",
    "        # i2 = 2 steps before (initial 2)\n",
    "        # i3 = 3 steps before (initial 3)\n",
    "        \n",
    "        ID = xyz1.finaldict[s][1]['id'][atom]\n",
    "        Type = xyz1.finaldict[s][1]['type'][atom]\n",
    "        \n",
    "        CNf0 = xyz1.finaldict[s][1]['c_cn'][atom]\n",
    "        CNi1 = xyz1.finaldict[step[s_index-1]][1]['c_cn'][atom]\n",
    "        CNi2 = xyz1.finaldict[step[s_index-2]][1]['c_cn'][atom]\n",
    "        \n",
    "        xf0 = xyz1.finaldict[s][1]['xu'][atom]\n",
    "        yf0 = xyz1.finaldict[s][1]['yu'][atom]\n",
    "        zf0 = xyz1.finaldict[s][1]['z'][atom]\n",
    "        rf0 = np.array([xf0, yf0, zf0])\n",
    "        nf0 = np.dot(rf0,a3)\n",
    "        xf0_prj = xf0 - yf0/tan    # Tilted projection of x\n",
    "        Nxf0 = math.floor(xf0_prj/Lx)    # x image number\n",
    "        Nyf0 = math.floor(yf0/Ly)    # y image number\n",
    "    \n",
    "        xf1 = xyz1.finaldict[step[s_index+1]][1]['xu'][atom]\n",
    "        yf1 = xyz1.finaldict[step[s_index+1]][1]['yu'][atom]\n",
    "        zf1 = xyz1.finaldict[step[s_index+1]][1]['z'][atom]\n",
    "        rf1 = np.array([xf1, yf1, zf1])\n",
    "        nf1 = np.dot(rf1,a3)\n",
    "        xf1_prj = xf1 - yf1/tan\n",
    "        Nxf1 = math.floor(xf1_prj/Lx)\n",
    "        Nyf1 = math.floor(yf1/Ly)\n",
    "        \n",
    "        xf2 = xyz1.finaldict[step[s_index+2]][1]['xu'][atom]\n",
    "        yf2 = xyz1.finaldict[step[s_index+2]][1]['yu'][atom]\n",
    "        zf2 = xyz1.finaldict[step[s_index+2]][1]['z'][atom]\n",
    "        rf2 = np.array([xf2, yf2, zf2])\n",
    "        nf2 = np.dot(rf2,a3)\n",
    "        xf2_prj = xf2 - yf2/tan\n",
    "        Nxf2 = math.floor(xf2_prj/Lx)\n",
    "        Nyf2 = math.floor(yf2/Ly)\n",
    "        \n",
    "        xi1 = xyz1.finaldict[step[s_index-1]][1]['xu'][atom]\n",
    "        yi1 = xyz1.finaldict[step[s_index-1]][1]['yu'][atom]\n",
    "        zi1 = xyz1.finaldict[step[s_index-1]][1]['z'][atom]\n",
    "        ri1 = np.array([xi1, yi1, zi1])\n",
    "        ni1 = np.dot(ri1,a3)\n",
    "        xi1_prj = xi1 - yi1/tan\n",
    "        Nxi1 = math.floor(xi1_prj/Lx)\n",
    "        Nyi1 = math.floor(yi1/Ly)\n",
    "        \n",
    "        xi2 = xyz1.finaldict[step[s_index-2]][1]['xu'][atom]\n",
    "        yi2 = xyz1.finaldict[step[s_index-2]][1]['yu'][atom]\n",
    "        zi2 = xyz1.finaldict[step[s_index-2]][1]['z'][atom]\n",
    "        ri2 = np.array([xi2, yi2, zi2])\n",
    "        ni2 = np.dot(ri2,a3)\n",
    "        xi2_prj = xi2 - yi2/tan\n",
    "        Nxi2 = math.floor(xi2_prj/Lx)\n",
    "        Nyi2 = math.floor(yi2/Ly)\n",
    "        \n",
    "        xi3 = xyz1.finaldict[step[s_index-3]][1]['xu'][atom]\n",
    "        yi3 = xyz1.finaldict[step[s_index-3]][1]['yu'][atom]\n",
    "        zi3 = xyz1.finaldict[step[s_index-3]][1]['z'][atom]\n",
    "        ri3 = np.array([xi3, yi3, zi3])\n",
    "        ni3 = np.dot(ri3,a3)\n",
    "        xi3_prj = xi3 - yi3/tan\n",
    "        Nxi3 = math.floor(xi3_prj/Lx)\n",
    "        Nyi3 = math.floor(yi3/Ly)\n",
    "    \n",
    "        dni1 = nf0 - ni1    # Change in normal from the past\n",
    "        dni2 = nf0 - ni2\n",
    "        dni3 = nf0 - ni3\n",
    "        \n",
    "        dnf1 = nf1 - nf0    # Change in normal in the future\n",
    "        dnf2 = nf2 - nf0\n",
    "        \n",
    "        # Test 1 step in the past\n",
    "        # Criterion: Change by more than 90% of interlayer distance\n",
    "        # History check: Must not revert back 1 or 2 steps in the past & in the future\n",
    "        if (abs(dni1) > 0.90*dn_avg) and (abs(dni2) > 0.90*dn_avg) and (abs(dnf1) < 0.10*dn_avg) and (abs(dnf2) < 0.10*dn_avg):\n",
    "            \n",
    "            ev_atom.append(atom)\n",
    "            dt = 1\n",
    "            \n",
    "            xf0 -= xy * Nyf0    # Wrap y first\n",
    "            yf0 -= Ly * Nyf0\n",
    "            \n",
    "            xi1 -= xy * Nyi1\n",
    "            yi1 -= Ly * Nyi1\n",
    "            \n",
    "            if (Nxf0 != 0) or (Nxi1 != 0):\n",
    "            \n",
    "                xf0 -= Lx * np.maximum(Nxf0,Nxi1)    # Shift initial & final x together toward the unit cell\n",
    "                rf0 = np.array([xf0, yf0, zf0])\n",
    "                nf0 = np.dot(rf0,a3)\n",
    "            \n",
    "                xi1 -= Lx * np.maximum(Nxf0,Nxi1)\n",
    "                ri1 = np.array([xi1, yi1, zi1])\n",
    "                ni1 = np.dot(ri1,a3)\n",
    "            \n",
    "            line = '%i\\t%i\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(s_index, dt, ID, Type, ni1, nf0, CNi1, CNf0)\n",
    "            out_ev.write(line)\n",
    "            \n",
    "            for aux in range(0,N_dump1):    # Detect any auxiliary atom that may move together\n",
    "                if aux not in ev_atom:    # Test only if not flagged before\n",
    "                    \n",
    "                    # Auxiliary atom labeled with \"i\" & \"f\" only\n",
    "                    \n",
    "                    ID = xyz1.finaldict[s][1]['id'][aux]\n",
    "                    Type = xyz1.finaldict[s][1]['type'][aux]\n",
    "                    \n",
    "                    CNf = xyz1.finaldict[s][1]['c_cn'][aux]\n",
    "                    CNi = xyz1.finaldict[step[s_index-1]][1]['c_cn'][aux]\n",
    "    \n",
    "                    xf = xyz1.finaldict[s][1]['xu'][aux]\n",
    "                    yf = xyz1.finaldict[s][1]['yu'][aux]\n",
    "                    zf = xyz1.finaldict[s][1]['z'][aux]\n",
    "                    rf = np.array([xf, yf, zf])\n",
    "                    nf = np.dot(rf,a3)\n",
    "                    xf_prj = xf - yf/tan\n",
    "                    Nxf = math.floor(xf_prj/Lx)\n",
    "                    Nyf = math.floor(yf/Ly)\n",
    "                    \n",
    "                    xi = xyz1.finaldict[step[s_index-1]][1]['xu'][aux]\n",
    "                    yi = xyz1.finaldict[step[s_index-1]][1]['yu'][aux]\n",
    "                    zi = xyz1.finaldict[step[s_index-1]][1]['z'][aux]\n",
    "                    ri = np.array([xi, yi, zi])\n",
    "                    ni = np.dot(ri,a3)\n",
    "                    xi_prj = xi - yi/tan\n",
    "                    Nxi = math.floor(xi_prj/Lx)\n",
    "                    Nyi = math.floor(yi/Ly)\n",
    "    \n",
    "                    dx = xf - xi\n",
    "                    dy = yf - yi\n",
    "                    dz = zf - zi                        \n",
    "                    dr = math.sqrt(dx**2 + dy**2 + dz**2)    # Change in position, rather than normal\n",
    "                    \n",
    "                    if dr > 0.80*dr_avg:    # Criterion: Change by more than 80% of the NN distance\n",
    "                        ev_atom.append(aux)\n",
    "                        \n",
    "                        # Auxiliary atom must be wrapped toward the active atom\n",
    "\n",
    "                        xf -= xy * Nyf    # Wrap y first\n",
    "                        yf -= Ly * Nyf\n",
    "                        \n",
    "                        xi -= xy * Nyi\n",
    "                        yi -= Ly * Nyi\n",
    "                        \n",
    "                        dxf = (xf - xf0)/Lx    # Extent of x deviation from the active atom\n",
    "                        dxi = (xi - xi1)/Lx\n",
    "                        \n",
    "                        if (abs(dxf) > 0.70) or (abs(dxi) > 0.70):    # Wrap if x deviation is more than 70% of the unit cell x dimension\n",
    "                                \n",
    "                            xf -= Lx * round(dxf)\n",
    "                            rf = np.array([xf, yf, zf])\n",
    "                            nf = np.dot(rf,a3)\n",
    "                            \n",
    "                            xi -= Lx * round(dxi)\n",
    "                            ri = np.array([xi, yi, zi])\n",
    "                            ni = np.dot(ri,a3)\n",
    "\n",
    "                        line = '%i\\t%i\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(s_index, dt, ID, Type, ni, nf, CNi, CNf)\n",
    "                        out_ev.write(line)\n",
    "        \n",
    "        # If no change 1 step in the past, test 2 steps in the past (change may be gradual)\n",
    "        elif (abs(dni2) > 0.90*dn_avg) and (abs(dni3) > 0.90*dn_avg) and (abs(dnf1) < 0.10*dn_avg) and (abs(dnf2) < 0.10*dn_avg): \n",
    "            \n",
    "            ev_atom.append(atom)\n",
    "            dt = 2\n",
    "            \n",
    "            xf0 -= xy * Nyf0\n",
    "            yf0 -= Ly * Nyf0\n",
    "            \n",
    "            xi2 -= xy * Nyi2\n",
    "            yi2 -= Ly * Nyi2\n",
    "    \n",
    "            if (Nxf0 != 0) or (Nxi2 != 0):\n",
    "            \n",
    "                xf0 -= Lx * np.maximum(Nxf0,Nxi2)\n",
    "                rf0 = np.array([xf0, yf0, zf0])\n",
    "                nf0 = np.dot(rf0,a3)\n",
    "            \n",
    "                xi2 -= Lx * np.maximum(Nxf0,Nxi2)\n",
    "                ri2 = np.array([xi2, yi2, zi2])\n",
    "                ni2 = np.dot(ri2,a3)\n",
    "            \n",
    "            line = '%i\\t%i\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(s_index, dt, ID, Type, ni2, nf0, CNi2, CNf0)\n",
    "            out_ev.write(line)\n",
    "    \n",
    "            for aux in range(0,N_dump1):\n",
    "                if aux not in ev_atom:\n",
    "                    \n",
    "                    ID = xyz1.finaldict[s][1]['id'][aux]\n",
    "                    Type = xyz1.finaldict[s][1]['type'][aux]\n",
    "                    \n",
    "                    CNf = xyz1.finaldict[s][1]['c_cn'][aux]\n",
    "                    CNi = xyz1.finaldict[step[s_index-1]][1]['c_cn'][aux]\n",
    "    \n",
    "                    xf = xyz1.finaldict[s][1]['xu'][aux]                            \n",
    "                    yf = xyz1.finaldict[s][1]['yu'][aux]\n",
    "                    zf = xyz1.finaldict[s][1]['z'][aux]\n",
    "                    rf = np.array([xf, yf, zf])\n",
    "                    nf = np.dot(rf,a3)\n",
    "                    xf_prj = xf - yf/tan\n",
    "                    Nxf = math.floor(xf_prj/Lx)\n",
    "                    Nyf = math.floor(yf/Ly)\n",
    "      \n",
    "                    xi = xyz1.finaldict[step[s_index-2]][1]['xu'][aux]\n",
    "                    yi = xyz1.finaldict[step[s_index-2]][1]['yu'][aux]\n",
    "                    zi = xyz1.finaldict[step[s_index-2]][1]['z'][aux]\n",
    "                    ri = np.array([xi, yi, zi])\n",
    "                    ni = np.dot(ri,a3)\n",
    "                    xi_prj = xi - yi/tan\n",
    "                    Nxi = math.floor(xi_prj/Lx)\n",
    "                    Nyi = math.floor(yi/Ly)\n",
    "    \n",
    "                    dx = xf - xi\n",
    "                    dy = yf - yi\n",
    "                    dz = zf - zi                        \n",
    "                    dr = math.sqrt(dx**2 + dy**2 + dz**2)\n",
    "                    \n",
    "                    if dr > 0.80*dr_avg:\n",
    "                        ev_atom.append(aux)\n",
    "                        \n",
    "                        xf -= xy * Nyf\n",
    "                        yf -= Ly * Nyf\n",
    "                        \n",
    "                        xi -= xy * Nyi\n",
    "                        yi -= Ly * Nyi\n",
    "                        \n",
    "                        dxf = (xf - xf0)/Lx\n",
    "                        dxi = (xi - xi2)/Lx\n",
    "                        \n",
    "                        if (abs(dxf) > 0.70) or (abs(dxi) > 0.70):\n",
    "\n",
    "                            xf -= Lx * round(dxf)\n",
    "                            rf = np.array([xf, yf, zf])\n",
    "                            nf = np.dot(rf,a3)\n",
    "\n",
    "                            xi -= Lx * round(dxi)\n",
    "                            ri = np.array([xi, yi, zi])\n",
    "                            ni = np.dot(ri,a3)\n",
    "                            \n",
    "                        line = '%i\\t%i\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(s_index, dt, ID, Type, ni, nf, CNi, CNf)\n",
    "                        out_ev.write(line)\n",
    "                            \n",
    "out_ev.close()\n",
    "\n",
    "t_end = time.perf_counter()\n",
    "print('Event detection done in %.2f' % (t_end - t_begin) + 's > event_flag.txt\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Event clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ev = open(file_ev,'r')    # Load list of time steps flagged as event points\n",
    "ev_lines = ev.readlines()\n",
    "ev_lines = [line.split() for line in ev_lines]\n",
    "\n",
    "for line_index, line in enumerate(ev_lines):\n",
    "    if line_index > 7:\n",
    "        line[0:4] = np.array(line[0:4]).astype(int)    # Time step, duration, atom ID, atom type\n",
    "        line[4:6] = np.array(line[4:6]).astype(float)    # Initial normal, final normal\n",
    "        line[6:8] = np.array(line[6:8]).astype(int)    # Initial coordination, final coordination\n",
    "\n",
    "ev.close()\n",
    "\n",
    "file_ls = pwd + '/event_list.txt'    # List of events\n",
    "out_ls = open(file_ls,'w')\n",
    "\n",
    "# List active atoms first, followed by event #, final time, and duration\n",
    "header1 = ' tf = final time (flagged step) \\n dt = event duration \\n ni = initial normal component \\n nf = final normal \\n CNi = initial coordination number \\n CNf = final coordination number \\n \\n'\n",
    "header2 = '\\t ID \\t type \\t ni \\t nf \\t CNi \\t CNf \\n event # \\t tf (ps) \\t dt (ps) \\n \\n'\n",
    "out_ls.write(header1)\n",
    "out_ls.write(header2)\n",
    "\n",
    "ev_no = 0\n",
    "ev_atom = []    # List of active atoms for each event\n",
    "ev_flag = []    # List of flags for each event\n",
    "\n",
    "for line_index, line in enumerate(ev_lines):\n",
    "\n",
    "    if line_index < 8:\n",
    "        continue\n",
    "    \n",
    "    elif line_index == 8:\n",
    "        \n",
    "        ev_flag.append(line)\n",
    "        ti = line[0] - line[1]    # ti = initial time\n",
    "        ID = line[2]\n",
    "        ev_atom.append(ID)\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        # If a step is within 3 steps of the previously flagged step, it is part of the same event\n",
    "        if line[0] - ev_lines[line_index-1][0] < 3:\n",
    "            \n",
    "            ev_flag.append(line)\n",
    "            ID = line[2]\n",
    "            \n",
    "            if ID not in ev_atom:\n",
    "                ev_atom.append(ID)\n",
    "\n",
    "            if line_index == len(ev_lines)-1:    # If at the last entry, finish clustering\n",
    "                \n",
    "                ev_no += 1\n",
    "                tf = ev_lines[line_index-1][0]\n",
    "\n",
    "                if tf == ti:\n",
    "                    dt = ev_lines[line_index-1][1]\n",
    "                else:\n",
    "                    dt = tf - ti\n",
    "\n",
    "                for atom in ev_atom:\n",
    "                    atom_flag = []    # List of flags for each atom\n",
    "                    \n",
    "                    for flag in ev_flag:\n",
    "                        ID = flag[2]\n",
    "                        \n",
    "                        if ID == atom:\n",
    "                            atom_flag.append(flag)\n",
    "                    \n",
    "                    for flag_index, flag in enumerate(atom_flag):\n",
    "                        \n",
    "                        ID = flag[2]\n",
    "                        Type = flag[3]\n",
    "                        ni = flag[4]\n",
    "                        nf = flag[5]\n",
    "                        CNi = flag[6]\n",
    "                        CNf = flag[7]\n",
    "                        dn = nf - ni\n",
    "                        \n",
    "                        if abs(dn) < 0.10*dn_avg:    # If normal does not change, print only if at the last flag\n",
    "                            \n",
    "                            if flag_index != len(atom_flag)-1:\n",
    "                                continue\n",
    "                            \n",
    "                            else:\n",
    "                                out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(ID, Type, ni, nf, CNi, CNf)\n",
    "                                out_ls.write(out)\n",
    "                        \n",
    "                        elif abs(dn) > 0.90*dn_avg:    # If normal changes, print that flag and move on to the next atom\n",
    "                            out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(ID, Type, ni, nf, CNi, CNf)\n",
    "                            out_ls.write(out)\n",
    "                            break\n",
    "\n",
    "                out = '%i\\t%i\\t%i\\n\\n'%(ev_no, tf, dt)    # Print the event number, time step, and duration\n",
    "                out_ls.write(out)\n",
    "                    \n",
    "        # If a step is after 3 steps or more from the previously flagged step, it is part of a different event, so finish clustering\n",
    "        else:\n",
    "\n",
    "            ev_no += 1\n",
    "            tf = ev_lines[line_index-1][0]\n",
    "\n",
    "            if tf == ti:\n",
    "                dt = ev_lines[line_index-1][1]\n",
    "            else:\n",
    "                dt = tf - ti\n",
    "\n",
    "            for atom in ev_atom:\n",
    "                atom_flag = []\n",
    "                \n",
    "                for flag in ev_flag:\n",
    "                    ID = flag[2]\n",
    "                    \n",
    "                    if ID == atom:\n",
    "                        atom_flag.append(flag)\n",
    "                \n",
    "                for flag_index, flag in enumerate(atom_flag):\n",
    "                    \n",
    "                    ID = flag[2]\n",
    "                    Type = flag[3]\n",
    "                    ni = flag[4]\n",
    "                    nf = flag[5]\n",
    "                    CNi = flag[6]\n",
    "                    CNf = flag[7]\n",
    "                    dn = nf - ni\n",
    "                    \n",
    "                    if abs(dn) < 0.10*dn_avg:\n",
    "                        \n",
    "                        if flag_index != len(atom_flag)-1:\n",
    "                            continue\n",
    "                        \n",
    "                        else:\n",
    "                            out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(ID, Type, ni, nf, CNi, CNf)\n",
    "                            out_ls.write(out)\n",
    "                    \n",
    "                    elif abs(dn) > 0.90*dn_avg:\n",
    "                        out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(ID, Type, ni, nf, CNi, CNf)\n",
    "                        out_ls.write(out)\n",
    "                        break\n",
    "\n",
    "            out = '%i\\t%i\\t%i\\n\\n'%(ev_no, tf, dt)\n",
    "            out_ls.write(out)\n",
    "            \n",
    "            ev_atom = []    # Reset for the next event\n",
    "            ev_flag = []\n",
    "            \n",
    "            ev_flag.append(line)\n",
    "            ti = line[0] - line[1]\n",
    "            ID = line[2]\n",
    "            \n",
    "            if ID not in ev_atom:\n",
    "                ev_atom.append(ID)\n",
    "\n",
    "out_ls.close()\n",
    "\n",
    "print('Event clustering done > event_list.txt\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LAMMPS NEB preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timestep = float(input('Enter the timestep used in the simulation, in ps : '))\n",
    "dmp_step = 0.1    # Dump frequency in ps\n",
    "d_step = int(dmp_step/timestep)    # Number of steps in a 0.1 ps interval\n",
    "\n",
    "# Header for LAMMPS DUMP file of initial trajectory\n",
    "l_step = 'ITEM: TIMESTEP'\n",
    "l_atom = 'ITEM: NUMBER OF ATOMS'\n",
    "l_box = 'ITEM: BOX BOUNDS xy xz yz pp pp ff'\n",
    "l_xyz = 'ITEM: ATOMS id type x y z'\n",
    "\n",
    "print('Preparing LAMMPS NEB images in ./neb/ev_***/ for each event...')\n",
    "\n",
    "ls = open(file_ls,'r')    # Load list of events\n",
    "ls_lines = ls.readlines()\n",
    "ls_lines = [line.split() for line in ls_lines]\n",
    "\n",
    "log = open(file_eq,'r')\n",
    "log_lines = log.readlines()\n",
    "\n",
    "for line_index, line in enumerate(ls_lines):\n",
    "    if line:\n",
    "        if (line_index > 9) and (len(line) == 3):\n",
    "        \n",
    "            ev_no = line[0]\n",
    "            tf = line[1]\n",
    "            dt = line[2]\n",
    "\n",
    "            tf_step = int(tf/timestep)    # Number of 0.1 ps steps\n",
    "            ti_step = tf_step - int(dt/timestep)\n",
    "\n",
    "            ev_dir = pwd + '/neb/ev_' + str(ev_no)    # Create NEB directory for each event\n",
    "            os.mkdir(ev_dir)\n",
    "\n",
    "            file_dmp = ev_dir + '/initial.dmp'    # LAMMPS DUMP file for visualizing initial trajectory\n",
    "            out_dmp = open(file_dmp,'w')\n",
    "\n",
    "            step = range(ti_step, tf_step+1, d_step)    # Event interval\n",
    "            for s_index, s in enumerate(step):\n",
    "\n",
    "                if s_index == 0:\n",
    "\n",
    "                    file_img0 = ev_dir + '/img0.dat'    # LAMMPS DATA file for initial structure (image 0)\n",
    "                    out_img0 = open(file_img0,'w')\n",
    "\n",
    "                    for l_index, l in enumerate(log_lines):\n",
    "\n",
    "                        if l_index <= head:\n",
    "                            out_img0.write(l)\n",
    "\n",
    "                    for atom in range(0,N_dump2):\n",
    "\n",
    "                        ID = xyz2.finaldict[s][1]['id'][atom]\n",
    "                        Type = xyz2.finaldict[s][1]['type'][atom]\n",
    "                        x = xyz2.finaldict[s][1]['x'][atom]\n",
    "                        y = xyz2.finaldict[s][1]['y'][atom]\n",
    "                        z = xyz2.finaldict[s][1]['z'][atom]\n",
    "\n",
    "                        out = '%i\\t%i\\t%f\\t%f\\t%f\\n'%(ID, Type, x, y, z)\n",
    "                        out_img0.write(out)\n",
    "\n",
    "                    out_img0.close()\n",
    "\n",
    "                header = '%s\\n%i\\n%s\\n%i\\n%s\\n%f %f %f\\n%f %f %f\\n%f %f %f\\n%s\\n'%(l_step, s, l_atom, N_dump2, l_box, xlo_ext, xhi_ext, xy, ylo, yhi, xz, zlo, zhi, yz, l_xyz)\n",
    "                out_dmp.write(header)\n",
    "\n",
    "                file_trj = ev_dir + '/coords.initial.' + str(s_index)    # Images 1 to N\n",
    "                out_trj = open(file_trj,'w')\n",
    "\n",
    "                header = '# ID \\t x \\t y \\t z \\n'\n",
    "                N_line = '%i\\n'%(N_dump2)\n",
    "                out_trj.write(header)\n",
    "                out_trj.write(N_line)\n",
    "\n",
    "                for atom in range(0,N_dump2):\n",
    "\n",
    "                    ID = xyz2.finaldict[s][1]['id'][atom]\n",
    "                    Type = xyz2.finaldict[s][1]['type'][atom]\n",
    "                    x = xyz2.finaldict[s][1]['x'][atom]\n",
    "                    y = xyz2.finaldict[s][1]['y'][atom]\n",
    "                    z = xyz2.finaldict[s][1]['z'][atom]\n",
    "\n",
    "                    out = '%i %i %f %f %f\\n'%(ID, Type, x, y, z)\n",
    "                    out_dmp.write(out)\n",
    "\n",
    "                    out = '%i\\t%f\\t%f\\t%f\\n'%(ID, x, y, z)\n",
    "                    out_trj.write(out)\n",
    "\n",
    "                out_trj.close()\n",
    "\n",
    "            out_dmp.close()\n",
    "\n",
    "ls.close()\n",
    "log.close()\n",
    "\n",
    "print('LAMMPS NEB preparation done > ./neb/ev_***/')\n",
    "print('Use initial.dmp to visualize the initial trajectory.')\n",
    "print('Use img0.dat as the initial structure.')\n",
    "print('Use coords.initial.** as the NEB images.\\n')\n",
    "\n",
    "print('Once the NEB calculations converge, NEB2DFT.py can be used to prepare DFT optimization.')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
